% Figure 3

% Axis options
axoptions = {
    'scaled ticks = false',...
    'y tick label style = {/pgf/number format/.cd, fixed, fixed zerofill, precision=3, set thousands separator={}}',...
    'x tick label style = {/pgf/number format/.cd, set thousands separator={}}',...
    'legend style = {font=\normalsize}'
}; 
  
% Load data
load ../data/DataforMatlab/IP.csv
load ../data/DataforMatlab/torn.csv

%% Percentile for plotting
run("../subroutines/sub_make_percentile.m")
%%Transform data to use boxchart functions
yearNum = zeros(10000,T);
for t = 1:T
    yearNum(:,t) = year(datetime(t+1973,1, 1));
end
yearAbrv = reshape(yearNum,[10000*T,1]);
 
% Vectorization
w_vectorize = reshape(w_vec_long,[10000*T,1]);

% Create table
tab2 = table(w_vectorize,yearAbrv);

% Calculation to find percentile considering household weight
hhw1974 = hhw_vec_L(:,:,1);
hhw2017 = hhw_vec_L(:,:,end);
logW1974 = log(Wfull(:,1));
logW2017 = log(Wfull(:,end));
hhw1974(isnan(logW1974)) = [];
hhw2017(isnan(logW2017)) = [];
logW1974(isnan(logW1974)) = [];
logW2017(isnan(logW2017)) = [];

% Rounding and maximum calculations for household weight
hhw1974 = round(hhw1974./10);
maxhhw1974 = max(hhw1974);
maxhhw2017 = max(hhw2017);

logW1974vec = NaN(size(logW1974,1),maxhhw1974);
logW2017vec = NaN(size(logW2017,1),maxhhw2017);

for i = 1:size(logW1974,1)
    logW1974vec(i,1:hhw1974(i)) = logW1974(i,1);
end
for i = 1:size(logW2017,1)
    logW2017vec(i,1:hhw2017(i)) = logW2017(i,1);
end

logW1974nan = reshape(logW1974vec,[1,size(logW1974,1)*maxhhw1974 ]);
logW1974nan(isnan(logW1974nan)) = [];

logW2017nan = reshape(logW2017vec,[1,size(logW2017,1)*maxhhw2017 ]);
logW2017nan(isnan(logW2017nan)) = [];

pct1974 = invprctile(logW1974nan,[log(min(U(:,T))) ; log(Uq(:,T)); log(max(U(:,T)))]);
pct2017 =invprctile(logW2017nan,[log(min(W(:,T))) ; log(Wq(:,T)); log(max(W(:,T)))]);
tt = 0:10:100;
vt = ones(size(tt))*pct2017(end);
vt2 = ones(size(tt))*pct2017(1);

%% Plot figures

figure('name','Figure 3','NumberTitle','off','Position',[100 100 1100 400])
subplot(1,2,2)
b = boxchart(tab2.yearAbrv,tab2.w_vectorize); hold on 

 xlim([1973,2018])
 b.JitterOutliers = 'on';
 b.MarkerStyle = '.';

upper = [NaN;log(Wmax)';NaN];
lower = [NaN;log(Wmin)';NaN];

plot(1973:2018,upper,'LineWidth',4,'Color',[0,0.5,0.9]); hold on %%[0,0.7,0.9]
plot(1973:2018,lower,'LineWidth', 4,'Color',[0,0.8,0.9]);  
legend('','Highest Utility','Lowest Utility','Location','northwest')
subtitle('Distribution of Total Expenditures over Time')
xlabel('Year')
ylabel('log Expenditure')
set(gca, 'FontName', 'Times')
% chenge font size==================
fontsize=11; 
h=gca; 
set(h,'fontsize',fontsize)

subplot(1,2,1)
plot(tt,tt,'--','Color','black','LineWidth', 2);hold on
plot(pct2017,pct1974,'Color',"#0072BD",'LineWidth', 4) 
plot(vt,tt,'--','Color','black','LineWidth', 1);
plot(vt2,tt,'--','Color','black','LineWidth', 1);hold off
yticks([0 20 40 60 80 100])
xticks([0 20 40 60 80 100])
 xlabel('2017 percentile')
 ylabel('Percentile of matched households in 1974')
subtitle('Matched Households in 1974 and 2017')
set(gca, 'FontName', 'Times')
% chenge font size==================
fontsize=11; 
h = gcf; %
set(gca,'fontsize',fontsize); %
%saveas(gcf,['../fig/Dists'],'epsc')
%exportgraphics(h,'../fig/Dists-eps-converted-to.pdf','BackgroundColor','none')




